CFDRC 


CFD Research Corporation 

3325 Triana Blvd. • Huntsville, Alabama 35805 • Tel: (205) 536-6576 • FAX: (205) 536-6590 


COMBUSTION CHAMBER ANALYSIS CODE 
Final Report 


ty 

A.J. Przekwas, Y.G. Lai, A. Krishnan, R.K. Avva, and M.G. Giridharan 


May 1993 


CFDRC Report: 4105/6 


for 

National Aeronautics and Space Administration 
George C. Marshall Space Flight Center 
Alabama 35812 


Contract Number: NAS8-37824 
NASA CORs: G. Wilhold and K. Tucker 


R&D Services and Software for Computational Fluid Dynamics (CFD) 


PREFACE 


Under the NASA MSFC contract NAS8-37824, CFDRC has developed an advanced 
CFD code: REFLEQS (REactive FLow EQuation Solver). REFLEQS has been 

designed for turbulent flow and heat transfer problems with and without chemical 
reactions. Particular attention has been given to the needs of predicting flow 
environment and efficiency of combustion in liquid rocket engines, using LOX/H 2 
or LOX/Hydrocarbon propellant systems. The results of this study are described in 
the following three volumes: 

1: Final Report; 

2: REFLEQS-Validation and Applications Manual; and 
3: REFLEQS - User's Manual. 

The final report includes descriptions of the mathematical basis of REFLEQS. 
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1 . 


INTRODUCTION 


1.1 Background 

Until recently, liquid rocket thrust chambers were designed using methods of 
characteristics and boundary layer approximations. In the mid 1980 s, new 
prediction methodologies, based on solutions of Navier-Stokes equations were 
proposed by Liang (1986), and Przekwas et al (1984, 1986). Computational Fluid 
Dynamics (CFD) found its way in practical rocket engine design due to several 
initiatives coordinated by NASA Marshall Space Flight Center and annually 
reported in the Proceedings of "Computational Fluid Dynamics Workshops". 

This report presents the results of the three-year effort (1989-1992), sponsored by 
NASA MSFC, to develop an advanced CFD code for predicting two-phase, reactive 
flows in liquid rocket thrust chambers. The starting point for this project was the 
existing REFLEQS code of CFDRC. 

1.2 Project Objectives and Approach 

The overall objective of this project was to develop a 3D CFD code that could be used 
as a tool in designing and predicting flow environment in LOX/H 2 and 
LOX/Hydrocarbon liquid propellant thrust chambers. The code should be capable of 
performing steady-state and transient analysis of single-phase and two-phase 
reactive flows on non-orthogonal grids. 

A building block approach was to be used in which the most advanced, 
methodology would be first implemented in the 2D code, tested and then extended 
to the 3D code. Code testing, verification /validation and documentation efforts 
were planned to be performed concurrently. 

Three major activities were planned for this project: development of the new 

generation REFLEQS code, code tests, and validation, and computational studies 
supporting STME engine development. 
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The starting version of the REFLEQS code utilized staggered grid arrangement, first 
order numerical schemes, interactive SIMPLEC algorithm and weak form of 
conservation equations (see Yang et at, 1992). The planned direction and selected 
code modifications for this project included the following: 

a. employ accurate differencing schemes (second-order accurate); 

b. incorporate efficient solution algorithms; 

c. incorporate options for Cartesian, cylindrical, and Body-Fitted 
Coordinate (BFC) systems; 

d. incorporate advanced turbulence models with boundary layer 
treatment; 

e. develop Eulerian-Lagrangian spray dynamics and combustion models, 
accounting for turbulence modulation effects for evaporating systems, 

f. incorporate available finite-rate reduced mechanism combustion 

models; 

g. incorporate turbulence-combustion interaction models; and 

h. include thermal radiation model. 

Early in the project, it was decided to incorporate state-of-the-art numerical 
methodology by implementing non-staggered grid arrangement, higher-order TVD- 
type numerical schemes, non-iterative PISO solution procedure, and strong form of 
conservation equations expressed in terms of Cartesian velocity components. 

Several new methodologies prior to this project had not been thoroughly tested, 
(e.g. PISO algorithm for compressible flows, Hautmann 4-step reduced chemistry for 
fuel rich conditions, etc.). As a result, fundamental development, testing, and 
validation efforts were directed to the development of new, robust fluid dynamic 
solvers. As a result, most code subroutines are completely reprogrammed. 

Parallel to the code development, active, continuous support was provided for 
analyzing flows and heat transfer problems during the STME development effort, 
e.g., NLS manifold study, STME nozzle cooling, vehicle base heating, etc. 
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1.3 Salient Fea tures of the Code 


This section outlines the salient features of the REFLEQS code, specifically designed 
for solving turbulent flow and heat transfer problems with and without chemical 
reaction. It consists of two fully compatible codes: REFLEQS-2D and REFLEQS-3D, 
and a common preprocessor, REFLEQSP. 

The codes are written in modular form, using FORTRAN 77 language. The 
provided pre-processor allows the user to specify the problem in a compact input file 
without modifying the main code. Several users can be supported with a single 
version of the code. Geometric flexibility implemented into the code allows the 
user to simulate flows in complex geometric configurations, multiple flow inlets 
and exits, and with internal solid objects. Selected capabilities, techniques, and 
physical models are summarized below. 

1.3.1 Numerical Techniques 

The following are numerical techniques of the REFLEQS code. 

a. Finite-volume approach of solving Favre-averaged Navier-Stokes 
equations. 

b. Cartesian, polar, and non-orthogonal Body-Fitted Coordinates (BFC). 

c. Colocated (non-staggered) grid. 

d. Strong conservative form of momentum equations with Cartesian 
components as dependent variables. 

e. Stationary or rotating coordinate systems. 

f. Pressure-based methodology for incompressible and compressible 
flows. 

g. Blockage technique for flows with internal solid objects. 

h. Four differencing schemes: upwind, central (with damping terms), 

MUSCL, and third-order Osher-Chakravarthy. 

i. Steady-state and transient (second-order time-accurate) analysis 
capability. 

j. Fully implicit and conservative formulation. 

k. Pressure-based solution algorithms, including SIMPLE, a variant of 
SIMPLEC, and non-iterative PISO algorithm for all flow speeds. 
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Symmetric Whole-Field linear equation solver and Conjugate 
Gradient Squared solver. 

Multi-component flows with heat and mass transfer. 
Eulerian-Lagrangian approach for two-phase flows. 


1.3.2 Physical Models 

The following are the physical models for the REFLEQS code. 

a. JANNAF Property Tables for gaseous species. 

b. Equation of State for arbitrary composition of ideal gas. 

c. Four turbulence models: 

1. Standard k-e model (Launder and Spalding, 1974); 

2. Extended k-e model (Y.S. Chen and Kim, 1987); 

3. Multiple Time Scale Model (Kim and C.P. Chen, 1988); and 

4. Low-Reynolds Number k-e Model (Chien, 1982); this provides 
the capability of resolving boundary layers. 

d. Standard and extended wall functions. 

e. Instantaneous reaction (diffusion controlled combustion) model. 

f. Equilibrium chemistry model (accelerated JANNAF method). 

g. One step finite rate chemistry models and two-step reduced 
mechanism models for 0 2 /H 2 and O 2 /Hydrocarbon propellants. 

h. Evaporating reactive spray dynamics model. 

i. Discrete ordinate thermal radiation model. 

1.3.3 Cndp Organization 

Figure 1-1 illustrates the overall flowchart of the code. A common preprocessor 
(developed under CFDRC's R&D funds) supports both 2D and 3D codes. It provides 
the flexibility of restarting problems, reading grids from external files, ( e.g ., EAGLE 
format or PLOT3D format), screening user's input file for simple errors and 
incorporates some degree of intelligence by analyzing input consistency. Results can 
be processed with NASA's PLOT3D or FAST graphical postprocessors, or more 
conveniently with CFDRC's CFD-VIEW scientific data visualization package. 
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Figure 1-1. REFLEQS Code Organization 


1.3.4 Code Testing 

The REFLEQS code has been checked out for a large number of problems by using 
the following four-step validation procedure: 

a. Check-out problems (with known solutions); 

b. Benchmark problems (with benchmark quality data); 

c. Validation problems (with detailed experiment data); and 

d. Field problems (with global performance data). 

Results of the computational studies can be found in the volume entitled, 
"REFLEQS - Validation and Applications Manual". 

1.4 Outline of the Rep ort 

Section 2 provides a general description of the governing partial differential 
equations embodied in the REFLEQS codes. Section 3 describes the physical models 
available. 
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The geometry and computational grid topologies are explained in Section 4 while 
the control volume discretization techniques are discussed in Section 5. Details 
including the similarities and differences of the solution algorithms, SIMPLE, 
SIMPLEC, and PISO, are described in Section 6. 

A separate section (Section 7) is dedicated to boundary conditions, providing 
detailed explanation for each boundary condition type available to the user. Section 
8 presents the overall solution procedure and linear equation solvers used in the 
code. Conclusions are presented in Section 9 and recommendations for further 
development and applications of REFLEQS are given in Section 10. 
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2 . 


BASIC GOVERNING EQUATIONS 


In this section, the basic governing equations for single-phase fluid flow are 
presented. These equations are derived from the conservation laws of mass, 
momentum, and energy which can be found in most fluid mechanics textbooks. 
Since a conservative finite-volume method is used in REFLEQS, all the governing 
equations are expressed in conservative forms. 

In Cartesian tensor form, the mass conservation and momentum conservation can 
be expressed as 


dp 

dt 




/ 


dp 


3t*. 


+ £( p ^)=-£: + ^ + <* 


(2.1) 


where w,- is the ith Cartesian component of the instantaneous velocity, p is the fluid 
density, p is the instantaneous static pressure, X\j is the viscous stress tensor, /; is the 
body force. For the Newtonian fluid, tjy can be related to the velocity through 


fdu { duj\ 


2 W 

3^ dx i. 


s a 


( 2 . 2 ) 


where ji is the fluid dynamic viscosity and is the Kronecker delta. 


The energy equation can take several forms and different forms are good for 
different classes of problems. Two forms are used in the REFLEQS code, static 
enthalpy; and stagnation enthalpy. The static enthalpy form of the energy equation 
can be expressed as: 


jp h 

dt 


d ( \ dtjj dp dp 

+ a^h V-W+dt^i^ 


T{ i dXj 


( 2 . 3 ) 
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Note that the above equation is not strictly conservative by its nature and is 
recommended for use for incompressible and low Mach number flows. The total 
enthalpy form of the energy equation, on the other hand, is fully conservative and 
is recommended for high speed compressible flows. The total enthalpy H, defined 
as H = h + v 2 j2, is governed by the following equation 


dp H d r rd W a .A 


tyj dp 


(2.4) 


In Equations (2.3) and (2.4), cjj is the j-component of the heat flux. By using Fourier 
Law, the heat flux is calculated by 


u dT 


(2.5) 


where k is the thermal conductivity. 

The species conservation equation is written as: 


apv f 

dt 





dYj\ 

dx } .) 


+ w 


( 2 . 6 ) 


where Y; is the mass fraction of species i, D is the mass diffusivity, and VV represents 
source/ sink terms due to chemical reactions. 

For turbulent flows, the diffusion terms in the above equations are replaced by the 
effective diffusion due to turbulence. Section 3 discusses the calculation of the 
effective diffusion (or eddy diffusivities) for various turbulence models. 
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3. 


PHYSICAL MODELS 


3.1 Turbulence Models 

In this section, the theoretical framework behind turbulence modeling is outlined 
followed by a detailed description of the various turbulence models available in the 
REFLEQS code. Section 3.1.1 introduces Favre time-averaging and the Boussinesq 
eddy viscosity concept. Section 3. 1.2-3. 1.7 describe the turbulence models 
implemented in the REFLEQS code. 

3.1.1 Favre Averaged Navier-Stokes Equations 

The basic governing fluid dynamics equations have already been introduced in 
Section 2. These equations are, in general, applicable to Newtonian fluid flow under 
steady or transient, incompressible or compressible, laminar, transitional or 
turbulent conditions. The nonlinearity of the Navier-Stokes equations, coupled 
with the complexity of the boundary conditions, makes it impossible to obtain 
analytical solutions for all but a limited number of flows of engineering interest. 
Hence one is forced to resort to approximate or numerical methods. Even though a 
wide variety of numerical techniques can be applied to solve the Navier-Stokes equa- 
tions for laminar flows. Direct Numerical Simulation (DNS) of turbulent flows is 
feasible only at very low Reynolds numbers. Turbulent flows are inherently 
unsteady and they contain a wide range of time and length scales, and resolution of 
these scales requires very short time steps and fine grids. The CPU and memory 
requirements for direct simulations are too large even for the fastest and largest 

present day computers. 


As most engineering applications only require time-mean quantities, the Navier- 
Stokes equations are usually averaged over time or ensemble of statistically 
equivalent flows to yield averaged equations. In the averaging process, a flow 
quantity <j> is decomposed in to mean and fluctuating parts. The following two types 
of averaging are generally used. 


Reynolds (or time) Averaging: 


r 

<{)=<}) + (}) where 



tydt 
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Favre (or density) Averaging: <J> = $ + ty” where $ = p<|>/p 

Note that overbar denotes Reynolds averaging while tilde denotes Favre averaging. 
T should be large compared the fluctuation time scale so that mean quantities are 
stationary over a number of samples. It also must be borne in mind that the mean 
quantities can vary in time on a scale much larger than T. 

Applying the Favre averaging procedure to Navier-Stokes equations (Equation 2.1), 
one obtains, the Favre-averaged Navier-Stokes (FANS) equations given below. (For 
detailed derivation, see Cebeci and Smith, 1974.) 




= 0 


dr dr dp 9 

aM + ^( pu ‘ u ;) = '^ + a^ 


it 


(&i 2 , 

k— + -r 2 - -T - — O,-. | 


du„ 


;L 


dX; dX; 3 dx 


m 


d£; b U i) 


(3.1) 


The FANS equations contain less information than the full NS equations, but have 

additional unknown terms - p u,u called the Reynolds stresses. These correlations 
between the fluctuating components arise in the averaging process, and need to be 
modeled to achieve closure of the FANS equations. 

All the turbulence models available in REFLEQS employ the generalized Boussinesq 

eddy viscosity concept in which the Reynolds stress - p iqw. is treated as a linear 
function of the mean strain rate 


du ; du ; 2 du 

1 , I 


p u i u i = ^\lk; + lki~3dxl 5{ i 




(3.2) 
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Here k is half the trace of the Reynolds stress tensor and fit is the turbulent eddy 
viscosity. Following the kinetic theory of gases, the eddy viscosity is modeled as the 
product of a velocity scale cj and a length scale 6. 


= C p (j (! (3.3) 

where C is a constant of proportionality. Various turbulence models differ in the 
way cj and p are estimated. In the following description of models, the overbar form 
jj. and p, and tilde for u, v, etc. are dropped for convenience. 

3.1.2 Baldwin-Lomax Model 

This belongs to the class of algebraic turbulence models because the velocity and 
length scales are obtained from algebraic relations. It is also commonly referred to as 
a mixing-length model because it employs Prandtl's mixing-length hypothesis in 
modeling length and velocity scales. 

Baldwin and Lomax (1978) developed this model primarily for wall-bounded flows. 
Like the mixing-length model of Cebeci and Smith (1974), it employs different 
expressions for /it in the inner and outer parts of the boundary layer. 

( M-t inner fo r V ~ V crossover 

| M-J outer for y ” y crossover 

In the inner layer, Prandtl's mixing-length model and the Van Driest's damping 
function are used to estimate the length scale, 

8 = Ky [l - exp (-y + / A + j] (3-5) 

where A+ is the Van Driest's damping constant. y+, the distance from the wall in 
wall units, is defined as 


(3.4) 


y + = yw T /v, m t = n /Vp 


(3.6) 
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In the preceding expression u T is commonly known as the friction velocity with t w 
being the shear stress at the wall. 

The velocity scale q is modeled as the product of f> and the root mean square 
vorticity. 



( 3 . 7 ) 


Using the preceding expressions, the eddy viscosity in the inner layer is obtained as 


V-tinner^P* 2 1 " 


( 3 . 8 ) 


The outer layer eddy viscosity is determined from the following expression: 


M'f outer ^ ^ cp P ^ wake ^ kleb (y) 

where K is the Clauser constant, is an additional constant, and 


( 3 . 9 ) 


F wake trnnl y max F max , C wk y max ^ 


u 2 ' 

u dif 




max 


The quantities ymax and F max are determined from the function 


( 3 . 10 ) 


( 3 . 11 ) 

The quantity F max is the maximum value of F(y) that occurs within the boundary 
layer and ymax is the value of y at which the maximum occurs. 


Fkieb(y) is the Klebanoff intermittency factor given by 
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F kleb ~ 



(Ckleby^ 

1 + 5.5 


^ V max ) 


-2 


( 3 . 12 ) 


The quantity Udif is the difference between the maximum and minimum total 
velocity in the boundary layer (i.e., at a fixed x station) 


Jmirt 


( 3 . 13 ) 


Utr </(“ 2 + 1,2 + + ^ + " ,2 )« 

The second term in Udif is zero for stationary walls. 

The values used for the constants appearing in the preceding expressions are 
A+ = 26 Cep = 1-6 Ckieb = 0-3 C w k = 0.25 k= 0.4 K = 0.0168 


3.1.3 Standard k-e Model 

Several versions of k-6 models are in use today, but the model employed in 
REFLEQS is based on Launder and Spalding (1974). This model employs two partial 
differential equations to estimate the velocity and length scales and hence 

commonly known as the two-equation model. These equations are the k- and e- 
equations which govern the transport of the turbulent kinetic energy (TKE) and its 
dissipation rate respectively. The modeled equations are 


d_ 

dt 


dx- 


, ,, „ 3 

> + dk 

(p“A)=p p -p £ + &; 

<r k dx L 


d_ 

dt 


dX: 


( nil r P Pe r P £2 , 3 

> + de 

(P u j E )- C z 1 Jt c e 2 k dXj 

°e dx L 


( 3 . 14 ) 


with the production P defined as 
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(du { dlij 2 du mA du i 2 d“r 
P - v ‘[to j + W i ~3dx m S ‘i)3x j ~3 k Sx„ 


(3.15) 


The square root of k is taken to be the velocity scale while the length scale is 
obtained from 


cW 2 


(3.16) 


The expression for eddy viscosity is 


V < = — 


(3.17) 


The five constants used in the model are: 


C^ = 0.09; C £j = 1.44 ; C z =132; a k =1.0; a z = 1.3 

Wall functions are used for this model. A description of wall functions is provided 
in Section 7.4. 


3.1.4 Extended k-e Model 

This model was developed by Chen and Kim (1987) by modifying the standard k-e 
model to make the dissipation rate more responsive to the mean strain rate than 
the standard k-e model. The k-equation retains the same form as that of the 
standard k-e model. In addition to the dissipation rate time scale, kje, a production 
range time scale, k/P is introduced in the e-equation as shown below. 


3 , , 3 ( \ r pPe r pe 2 , r PP 2 . 3 3e 

*b £ ) + ^[ pu i E r c e—- c *— +c z 3 — + ^ pra^j 


(3.18) 
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The net effect of this formulation is to enhance the development of e when the 
mean strain is strong and to suppress the generation of e when the mean strain is 
weak. The model constants are: 


C„ = 0.09 ; c e =2.25 


C 


. =1.9; C F =0.25; o k = 0.75; a t = 1.15 

-2 fc 3 * 


3.1.5 Multiple Time-Sca le Model 

Most turbulence models, including the standard k-e model, assume one time scale 
for both the production and the dissipation rates of turbulence. However, experi- 
ments and the Direct Numerical Simulation (DNS) of turbulence have shown that 
most of the turbulence production occurs at large scales (energy carrying eddies) and 
is cascaded to smaller scales (dissipative eddies) where most of the dissipation takes 
place. Several authors (e.g., Hanjalic et ah, 1980; Fabris et ah, 1981) have suggested 
partitioning of the energy spectrum into production range k p and dissipation range 
k t/ the sum of the two quantities being k. Kim and Chen (1988) have developed a 
multiple time-scale turbulence model based on variable partitioning of turbulent 
kinetic energy. The partitioning is determined as a part of the solution and depends 
on local turbulence intensity, production, energy transfer and dissipation rate. The 
partition is moved into higher wave numbers when production is high, and to low 
wave numbers when production vanishes. The model uses a transport equation for 
each k p , k t , e p and e t , however these equations differ from each other only in the 

source terms. 


lM + s;( p “;*p) =pP ‘ p vs; 

lM + 4-( p ^) =p v pe < + s; 


f p + dk p \ 

* i) 
f\i + [L t dk t \ 

[ °kt dx ]j 


P P 2 


pP% 


dt ( p£ p) + dx Ce Pi kp +( ~ £ P2 kp kp 


P e o 


dXj 


l ^ fyj 
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(3.19) 


3. , 3 , . r PeP 2 . c gfgpfi c + 

jN^W'^TT Cb . kt “• kt 3*, ( 00 , dXj) 

The eddy viscosity is defined as 


U( = C M fc 2 /^ 


( 3 . 20 ) 


where k, the total TKE, is calculated as 


k = k p + k t 


( 3 . 21 ) 


The model constants are: 


Cfi = 0.09; (Jkp — <Jkt — 0.75; — <Jg 1.15 

Cepi =0.21; Cep! =1.24; C £p3 = 1.84;C e n =0.29; C & 2 =128; ^=1.66 


3.1.6 Low Reynolds Number Model 

All three models described above, viz. the Standard k-e model, the extended k-e 
model, and the multiple time scale model, are of high-Reynolds number form and 
therefore require the used of wall functions. However, the commonly used wall 
functions may not be accurate in flows with large separation, suction, blowing, heat 
transfer, relaminarization, etc. This difficulty associated with wall functions can be 

reduced by the use of low-Reynolds number k-e models which permit integration of 
momentum and k-e equations all the way to the wall. The k-e equations are 
modified as shown below to include the effect of molecular viscosity in the near 

wall regions. 
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s ( p*) + e-(p «,•*) = £; 


d 

H + lh dk 


& 

J*1 


d f , a 

(pe) + 


3f 




1- 9 

de 

i “ dXj 



p(P-e-D) 


+ C e ,/,^-C £ i / 2 ^ + E 


- , P * 2 

Pi = c p/^T 


( 3 . 22 ) 


Several versions low-Reynolds number k-e models are available today. After a 
careful analysis of the comparative studies done by Patel et al. (1985), Brankovic and 
Stowers (1988), Awa et al (1990), the low-Re model of Chien (1982) was selected for 
implementation in REFLEQS. The model parameters appearing in the preceding 
equations are: 


= 0.09 ; C £i = 135; C t =lJ8 ; < 7 * = 1.0 ; o e = 1.3 
f^l-exp(-0.011Sy*),f 2 = 1.0; f 2 = l -022exp[{RJ6f] 
D = 2 vk/y 2 ; E = -2v{e/y 2 ) exp (-0-5 y + ) 


3.2 Combustion Models 

The combustion models implemented in the REFLEQS code include: 
instantaneous, equilibrium, and finite rate chemistry models. Previously, the finite- 
rate model in the REFLEQS code was restricted to one-step or multi-step irreversible 
reactions such that a particular fuel species may appear as a reactant in one reaction 
step only. In this study, efforts were made to extend the chemistry modeling 
capability of the REFLEQS code to include dependent, multi-step reactions in which 
the same fuel (or oxidizer) species can exist in more than one step. The cases of H 2 - 
Air and Hydrocarbon- Air combustions are considered. 

3.2.1 2-Step H r O? Reaction Model 

The combustion model for H 2 -O 2 reaction is based on the 2-step model of Rogers 
and Chinitz (1983). 
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H 2 + 0 2 = 20H 
H 2 +20H = 2H z O 


(3.24a) 

(3.24b) 


Further simplifications are made by assuming that the first step. Equation (3.24a), is 
in equilibrium, and the second step. Equation (3.24b), is irreversible. The 
implementation of this reaction mechanism into the REFLEQS code is discussed 

below. 

Basic Equations: The basic approach of REFLEQS in solving chemical species, is to 
solve the mixture fractions and the // fuel ,/ species of a finite-rate reaction first, and 
then derive the mass fractions (or concentrations) of all species according to element 
conservation. For this 2-step H 2 -O 2 reaction of mixed equilibrium chemistry and 
finite-rate chemistry, the basic equations are the tranport equations of conserved 
mixture fractions 


F (fi) = 0 , i = 2, 2, ..., NCMPS - 1 (3-25) 

where NCMPS is the number of mixture compositions. The transport equation of 
fuel species [H 2 ] in Equation (3.24b) is given as: 

fN-VpNIo ' 1 ! 2 (326) 

The extra relations include equilibrium Equation of (3.24a), i.e., 


\ OHf 

K °h-[H 2 ][0 2 ] 


(3.27) 


and element conservation 


[0H] + 2[H 2 ] + 2[H 2 0] = [H] 


(3.28a) 
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(3.28b) 


[0H] + 2[0 2 ] + [H 2 0\ = [0] 

In the above Equations (3.25) through (3.28a,b), F( ) denotes the fluid transport 
operator and,/; is the mixture fraction to be solved for NCMPS-1 compositions (with 

NCMPS-l 

fNCMPS =1 ~ " /»). The equilibrium constant K,oh for Equation (3.24a) is 

expressed as 


Kqh ~ ex P 


, AG - 2G 0 h" ( g h 2 + G Oj) 


(3.29) 


where G is the Gibbs free energy at T. 

Solution Procedure: Once the mixture fraction /; is solved, the mole values of 

elements [H], and [O] can be derived as 


[H]=[0Hf , + 2[H 2 f' + [H 2 0f' 

(3.30a) 

[0] = [OH]™ + 2 [O 2 ] <0) +[H 2 of ' 

(3.30b) 


where [ ](© is the un-reacted species concentration from mixture fraction /,-. 

With [ H 2 ] solved from Equation (3.26), the species [Qzl, [OH], [H 2 O] are to be obtained 
from Equation (3.30b) plus 


[0H] + 2[H 2 0] = [H]-2[H 2 ] 


(3.31) 


[OH] 2 = K oh [h 2 ] [o 2 

After algebraic manipulations, the Equations (3.30b), 
different cases of [H 2 ] values as: 


(3.32) 

(3.31), and (3.32) are solved in 
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[H 2 ] > 0, then 


(i) 


[OH] = 1(7? 

+ 4bRj -b) 

(3.33a) 

[oj- my a 


(3.33b) 

(ii) [H 2 ]= 0 , then 



[OH] = 0 
[ 0 2 ] = J/4Rj 


(3.34a) 

(3.34b) 

and in both cases 

[H 2 0] = [0]-2[0 2 ]-[0H] 

(3.35) 

In the above expressions. 

b = l/4a = l/4K OH [H 2 ] 

(3.36) 


Rj = max (o, 2[0] - [H] + 2 ^H 2 jJ 

(3.37) 

3.2.2 2-§tep CHj. - O 2 Reaction Model 

The 2-step CH 4 - 0 2 reaction model of Westbrook and Dryer (1981) is 
simplified by assuming only forward reaction is present. 

employed and 


£ 

CH 4 + ^0 2 ^C0 + 2H 2 0 

(3.38a) 


2 

—0 2 + C0 CO 2 

(3.38b) 

where the reaction rates of Equations (3.38a and 3.38b) are given by Westbrook and 
Dryer (1981), (in cm-g-s unit) respectively as 


20 


4105/6 


(3.39a) 


q = 6.7 x 10 12 e~RT~~ [CH 4 ] 02 [O 2 ] 2 ' 3 

d* = 10 24 ' 6 g TT [CO] [H 2 0]° 5 [0 2 f 25 (3.39b) 

An algorithm similar to that for the 2-step H 2 -O 2 mechanism is used. The transport 
equations to be solved are Equation (3.25) for mixture fraction fi and those for "fuel" 
species as 


FflCHj)^ (3 - 40) 

F J0 2 J = -1 - 056^ (3.41) 

Other species can be derived from the element conservation relations. 

[CH 4 ] + [CO] + [C0 2 ] = [C] (3.42a) 

4[cH 4 ] + 2[h 2 0] = [H ] (3.42b) 

2 [0 2 ] + [CO] + 2 [C0 2 ] + [h 2 0] = [O] (3.42c) 


With elements [C], [H], and [O] obtained from mixture fraction fi, and [ CH 4 ], [0 2 ] 
solved from Equations (3.39), (3.41), and (3.42) yields 


[H 2 0} = t{H]-2[cH t \ 

(3.43) 

[C0 2 ] = [0] - [C]+ [CH 4 ] - 2 [0 2 ] + [HjO] 

(3.44) 

[co]=[c]-[ch 4 ]-[co 2 ] 

(3.45) 
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3.2.3 Chemical Eq uilibrium Model 

The equilibrium model proposed by Meintjes and Morgan (1989) is described below. 
This method is based on calculating the element potentials or the element activities. 

The chemical potential of any species can be written as: 


p r p° j + RTln 



(3.46) 


where Uj = chemical potential, p° = standard state Gibbs function, p 0 - atm • pr and 
Pj = partial pressure of species " j ". 



(3.47) 


where rij/n is the mole fraction of species p is the total pressure of mixture and 
p 0 is atm pressure. 

The chemical potential "pf' can be written as: 

I 

Pj = L a {j M- f ( 3 . 48 ) 

1 = 2 

where pi is the chemical potential of atomic element i and fly is the number of moles 
of element "i" per mole of species Substituting into Equation (3.46), 


1 * _ V) . fap) 

RT Xfl;yM-f RT + n [n p 0 j 


(3.49) 


From the equation of state. 


p _RT 
n ~ V 


(3.50) 
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Therefore, 


RT] 




Hi 


1 = 2 


Vo 


V RT 


Defining molar concentration, cj = tijfv and 


Rewriting Equation (3.51), 


. V°j , (RT\ i i 
lnc i~ RT "(pj + 


where Cj = wy/v and A,- = fli/RT. 
Equation (3.52) can be written as: 


In Cj = In Rj + 2 &ij In Z ( - 

i= i 


where 


jy PO 

and Z,- = (an element variable) 

Equation (3.53) is simplified as: 

c r R / n(z I p/-i,.../ 

t=2 



(3.51) 


(3.52) 


(3.53) 


(3.54) 


(3.55) 
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Since atoms of each element are conserved, 

i 

Z a i j c j = b i/ i = l,...I (3.56) 

/=2 

where bi is the total number of moles of element "i" in all species. 

Equations (3.55) and (3.56) form the basis of this formulation. Substituting Equation 
(3.55) into (3.56): 


iL j R j n{z m )A-b i =o,i=i,...i 

j=l{ m=l ) 


(3.57) 


Thus there are "i" non-linear algebraic equations for the "I" unknown values of Z,-. 
A good guess is used to obtain the initial values of Zf and the subsequent values are 
obtained by solving this system. The values of cj are calculated from Equation (3.55). 
The temperature is calculated from “cf and then the new values of Rj are calculated 
and again the system is solved for new values of Z; and so on until convergence is 

achieved. 


3.3 Turbulence-Combustion Interaction 

In many combustion systems, reactants enter in separate streams, with fuel and 
oxidizer non-premixed. The resulting combustion process of diffusion flames, in 
the presence of flow fluctuation due to turbulence, is very complicated. Generally, 
the "fast-chemistry" assumption is made to simplify the problem of turbulence- 
combustion interaction. Under this assumption, the chemical reactions are 
considered to have very high production rates, so that the reactions are completed as 
soon as the reactants are mixed. In other words, the reaction time is assumed to be 
negligibly short in comparison with the turbulent mixing time. In the limit of fast 
chemistry, the instantaneous species concentrations and temperature (with given 
enthalpy) are functions only of conserved scalar concentrations (s.g., mixture 
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fractions) at that instant. The averaged thermodynamic variable can be obtained 
from the knowledge of statistics of those scalars, which can be expressed by the 
probability-density function (pdf). A further approximation is made by assuming a 
form of the pdf that can be evaluated from a few moments of the conserved scalar. 


In this study, the conserved-scalar approach based on the fast-chemistry assumption 
is adopted to account for the effects of turbulence-combustion interaction. The pdf 
of the mixture fraction is assumed to be a tophat function. 

3.3.1 Conserved Scalar Equations 

For the problem of diffusion controlled combustion, the effect of mass diffusivity 
differences among different species are negligible in most turbulent flows at 
moderate to high Reynolds number (see Kuo, 1986). Therefore, the simplest and 
best choice for the conserved scalar varaible is the mixture fraction. 


To determine the pdf, which describes the statistics of mixture fraction, the mean 
value and variance of the mixture fraction are needed. These can be solved from 
the Favre-averaged transport equations suggested by Jones and Whitelaw (1982), as 


df dptdf) 

pU jd^- = ^\T t Wj (3-58) 


- ~ dg d fM-t dg\ 
pu id Xj dXj \c f dxjj 


+ 2 


hfdf 


?\2 


r Pe~ 

~ C DT 8 


(3.59) 


where f and g (g = f" 2 ) are the density- weighted mean and variances of mixture 
fraction/, respectively, and Cq has the value of 2. 

3.3.2 Top-Hat Pdf Model 

The "top-hat" pdf with possible delta function at 0 and 1 is assumed as 
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(3.60) 


f P o ' f=° 

pif) = lc, 0<a </< b<l 

\Pi'f=l 


The specific expression of p(f) ( i.e the values of parameters, P 0 , Pj, a, b, c ) can be 
determined by the following moments. 


[pWM 

(3.61a) 

[ fp(f) d f=f 

(3.61b) 


(3.61c) 

which leads to the explicit form for Equation (3.61) as 


P 0 + P 1 +c(b-a) = l 

(3.62a) 

P 1 + ^-c(b 2 -fl 2 ) = / 

(3.62b) 

P 1 +^c{b 3 -a 3 ) = g + f 2 

(3.62c) 


The expression of p(f) for the mean value f < 0.5 can be derived for three different 
cases from Equations (3.62a-c) as 
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Owing to the symmetrical properties of and anci 

f = 0.5 for Je(0,i) , for f >05, P(f) can be expressed simply by obtaining the 
parameters as 


a\j= 2-6 1( 2 -/) (3.66a) 

b\j= 2— |(2 _ /) (3.66b) 

c \f = C \{i-J) (3.66c) 

P o\f~ P l |( 2 _/) (3.66d) 

^i|/ = %_/) (3.66e) 


where []|/ denotes the parameter for f >05, and [ ] |(i - /) is the parameters given 
in Equations (3.63-3.65) with j substituted for (l - f) . 

3.3.3 Mean Values of Thermodynamic Variables 

Based on the fast-chemistry assumption, the thermodynamic variables such as 
species mass fractions (or concentrations), temperature, and density, are functions 
only of the mixture fraction. The mean values of these variables, which are of 
interest, are obtained by averaging the instantaneous values at the given mixture 

fraction, /, with the pdf P(f) over the sample domain fe[0 ,2} i.e. 


Yi-fYiWWif 

Jo 


(3.67a) 




(3.67b) 
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(3.67c) 



in which Y; is the species mass fraction, T, p are temperature and density, 
respectively. 

The integration of Equations (3.67a-c) are approximated by numerical quadrature 
with limited discrete samples. 

3.4 Sprav Model 

The droplet equations of mass, momentum and energy are solved in a Lagrangian 
frame of reference moving with the droplets. These solutions are used to calculate 
the source/ sink terms for the corresponding gas phase equations. 

The equation of motion for the droplet is written as (Crowe et al. (1977)): 

m d % = C D p(U-v)|U-v| — -V d Vp + m d g ( 3 . 68 ) 

where m d is the mass of the droplet and v = ui + v\ + zvk its velocity vector, u, v 
and w are the Cartesian velocity components. C D is the drag coefficient and p, U 
and p are the density, velocity and pressure of the gas, respectively. A d is the droplet 
surface area and V d is the droplet volume. For a spherical droplet, A d = mV and V d 

= nd3/6 where d is the droplet diameter, g is the gravity vector. Equation (3.68) 
accounts for the acceleration /deceleration of the droplet due to combined effects of 
drag with gas, local pressure gradients and body forces such as gravity. 

The drag coefficient, C D , for the droplet is based on the local Reynolds number 
evaluated as: 


\ 
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(3.69) 


Re = Pj£l Ai 

P 

where p is the dynamic viscosity of the gas. The following correlations have been 
found to be valid for different ranges of Reynolds numbers (Crowe et ul, 1977). 

Cn = — for Re < 1 

Re 

C D = -2![l + 0.15Re°- 687 ] for 1 < Re < 10 3 

Re (3.70) 

C D = 0.44 for Re > 10 3 


Substituting for Cd into Equation (3.68), 


dv 

It 


ps 2 


(U 


, V P 

v)-^ + g 


(3.71) 


where pd is the droplet liquid density and /= CoRe/24. The droplet locations are 
determined from the velocities as: 


dK = v (3.72) 

dt 

where R is the position vector of the droplet. 

The mass conservation equation for the droplet is given by: 

^±=-Sh(pD)7d(x v -x 00 ) (3.73) 

where Sh is the Sherwood number and D is the mass diffusion coefficient for the 
gas. x v and Xoo are the vapor mass fractions at the droplet surface and the free stream, 
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na i Eli b 


respectively. In this model, x v is calculated from the liquid saturation pressure. The 
Sherwood number is obtained from the following expression (Crowe et ul, 1977). 

Sh = 2 + 0.6Re°- 5 S c 0333 (3,74) 

where S c is the Schmidt number. Since m d = p d 7ui 3 /6, the mass equation for the 
droplet is rewritten in terms of its diameter: 


M_. 2Sh(pD)^i 
dt v ' Pdd 

(3.75) 

The energy equation for the droplet is written as: 

+ 0.76) 

where q is the sensible heat transferred to the droplet and T d is the droplet 
temperature. L is the latent heat of vaporization for the droplet fluid and C d is the 
specific heat of the droplet, q is calculated as: 

cj = Nu jckd (Tg - T d ) (3.77) 

where k and T g are the thermal conductivity and temperature of the gas, 
respectively. The Nusselt number, Nu, is obtained from the following correlation: 

Nu = 2 + 0.6Re°- 5 Pr^- 333 (3 ' 78) 

where Pr is the Prandtl number. Substituting the above expressions into Equation 
(3.76), the energy equation is rewritten as: 


dTj _ T g ' T _1 9k 

dt o d 

where 


(3.79) 
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(3.80) 


LSh(pD)(x P -x 00 ) 
® L ~ Nuk 


0 = 


p^c s 

6Nuk 


(3.81) 


If the droplet temperature is below its boiling point. Equation (3.75) is used as the 
mass conservation equation. Once the droplet attains its boiling point, the mass 
conservation equation is obtained by setting the left hand side of Equation (3.76) to 
zero, i.e., the droplet temperature cannot exceed its boiling point. Therefore, 


drEd = . I (3.82) 

dt L 

This equation implies that all the heat transferred from the gas is used to vaporize 
the droplet. 

The solutions to the droplet equations provide sources/ sink terms to the gas phase 
conservation equations. The sources/ sinks are assigned to the particular cell (in the 
Eulerian gas phase) in which the droplet is located. The calculation of the 
source/sinks is as follows. 

The mass of the droplet is continuously monitored along its trajectory. The 
source/ sink term for a cell is given as: 


Am = np d 7t X 

Z 



(3.83) 


the summation takes place over all the droplets located in the given cell. The 
subscripts "in" and "out" refer to cell inlet and outlet conditions, n is the number 
density of droplets in a given parcel. Similarly the momentum source/sink term is 
evaluated as: 
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(3.84) 


AM = np d K X 


ML'MiL 


This includes the effects of frictional interactions with the gas. The momentum 
source term is a vector. The energy source/ sink term is given as: 


AE = Amh v -n'£i (3.85) 

i 

where h v is the enthalpy of the vapor at the vaporization temperature and q. is the 
heat transfer from the gas to the droplet (calculated according to Equation (3.77)). 

3.5 Radiation Model 

A large number of numerical techniques are available for solving the equations 
governing the transfer of thermal radiation. The most widely used among them are 
listed below. 


a. Hottel's Zone method 

b. Flux method 

c P-N approximation method 

d. Monte-Carlo method 

e. Discrete-transfer method 

f. Discrete-ordinate method 

Hottel's zone method was the first general numerical procedure developed to 
handle the radiative transfer equation. It has been widely used to estimate the 
radiation transfer in the absence of detailed knowledge of the flow and reaction. 
This method involves the use of exchange factors which needs to be worked out in 
advance for complex geometries. The disadvantage of this method is that it is very 
uneconomical and not suited for complex configurations with strongly participating 
media. The flux methods involve differential approximation to the radiative 
transfer equation. This method is computationally efficient and for this reason it 
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has been widely used in general combustion problems. However, this method 
suffers from inadequate accuracy and the difficulty to extend to complex geometries. 
Compared to this method, the P-N approximation method is more accurate. 

Among all the existing numerical treatments for radiative heat transfer, the Monte- 
Carlo method is the most accurate method. In this method, the randomly chosen 
gjigrgy releases are tracked through the combustion chamber for their lifetimes. 
Unfortunately, this method is computationally expensive. The discrete-transfer 
method, on the other hand, is economical and can be applied in a straightforward 
manner to complex geometries. This method, however, cannot treat the scattering 
by the medium containing droplets and soot particles and hence is not suitable for 
combustor applications. 

After comprehensive review of all the methods, it has been found that the discrete- 
ordinate is well suited for radiation prediction in rocket combustion chambers. 
There are several advantages of this method. Firstly, the discrete-ordinate method 
can be easily extended for polar and complex geometries. Secondly, the evaluation 
of the in-scattering term is relatively simple. Also, the boundary conditions can be 
imposed more accurately at inlets and exits. Finally, this method is accurate and can 
be easily coupled with CFD codes based on the control volume approach. The 
governing equations and the solution procedure for this method are described 
below. 

3.5.1 Governing Equations 

The integro-differential radiative heat transfer equation for an emitting-absorbing 
and scattering gray medium can be written as (Fiveland, 1984) 

(Q-V)/(r,Q) — (k + oJffr.QJ + WjM + i f. I(r,fll 4 (a'-t 0 )ja (3.86) 


where Q is the direction of propagation of the radiation beam, I is the radiation 
intensity which is a function of both position (v) and direction (i2), k and o' are the 
absorption and scattering coefficients respectively. is the intensity of black body 
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radiation at the temperature of the medium, (j) is the phase function of the energy 
transfer from the incoming Q' direction to the outgoing direction Q. The term on 
the left hand side represents the gradient of the intensity in the specified direction 12. 
The three terms on the right hand side represent the changes in intensity due to 
absorption and out-scattering, emission and in-scattering respectively. 

The boundary condition for solving the above equation may be written as 

I(r,£2) = el b (r) + | f , \n-ri\l[r,a)dri (3.87) 

J n ■ £1 < 0 


where I is the intensity of radiant energy leaving a surface at a boundary location, e 
is the surface emissivity, p is the surface reflectivity, and n is the unit normal vector 
at the boundary location. 

In the discrete-ordinate method. Equations (3.86) and (3.87) are replaced by a discrete 
set of equations for a finite number of ordinate directions. The integral terms on the 
right hand side of the above equations are approximated by summation over each 
ordinate. The discrete-ordinate equations may then be written as. 


dL 


dL 


m, 3^ + - -I (*• + °) L + + m I. »„■' I m - ; m = 1. . .. M 


% 


(3.88) 


m 


For cylindrical polar geometries, the discrete-ordinate equation may be written as. 


5m d( r jj I _ j_ a (n Jm[ d J^ = _ [k+a)I +R+^-Ia) '(|) - = 

~ dr r 34> M ' m dz 1 J m b 4TI m Vmm m 

771 


(3.89) 


where the second term on the left hand side represents the angular intensity 
redistribution. It accounts for the change in the direction cosines as a beam travels 

in a straight line. 


35 


4105/6 


The different boundary conditions required to solve this equation are given as. 


Wall Boundary 



I m = d b + Yl^ W m\^m\ I m 
m 

(3.90) 

Symmetry Boundary 


^ m ' 

(3.91) 

Inlet and Exit Boundary 


7=0 

m 

(3.92) 


In the above equations, m and m' denote the outgoing and incoming directions, 

respectively. For a particular direction, denoted by m, the values fi and £ are the 
direction cosines along the x and y directions. These equations represent m coupled 
partial differential equations for the m intensities, I m . In the present code, only the 
S 4 approximation, i.e., 12 ordinate directions are considered. Table 3-1 lists the 
values of the direction cosines and the associated weights. 


Table 3-1. 

Direction Cosines and Their Associated Weights for the S 4 Approximation 


Direction 

Number 


Direction Cosines 


Weight 

w m 

1 

-0.908248 

-0.295876 

.295876 


2 

-0.908248 

0.295876 

.295876 

jt/3 

3 

-0.295876 

-0.908248 

.295876 

71/3 

4 

-0.295876 

-0.295876 

.908248 

71/3 

5 

-0.295876 

0.295876 

.908248 

71/3 

6 

-0.295876 

0.908248 

.295876 

71/3 

7 

0.295876 

-0.908248 

.908248 

71/3 

8 

0.295876 

-0.295876 

.908248 

tc/3 

9 

0.295876 

0.295876 

.908248 

71/3 

10 

0.295876 

0.908248 

.295876 

71/3 

11 

0.908248 

-0.295876 

.295876 

71/3 

12 

0.908248 

0.295876 

.295876 

71/3 

4105/6 13-1 


36 


4105/6 



Equation (3.88) can be cast in a fully conservative form for a general BFC coordinate 
system as 


^7 [4 m imJFlx + ^mLJ F ly] K l m l F 2x + 7 m J F 2y] 


= / 


— (fc + 8) I m + k 1^ + 2 


(3.93) 


where x' and y' are the local grid coordinate directions, / is the Jacobian of the 
coordinate transformation and 


F 

lx dx 
F 

F iy~dy 


_ V 

2* 3x 


■ = 2L 

2y 


(3.94) 


Equation (3.89) can be integrated over the control volume, shown in Figure 3-1, to 
obtain 


M-m (A e Im t A w I m Jj + £> m I m ^ A s ^s) 


[fm + lf2h,mt-ll2 I'm -7/2 


W 


m 


— - (fc + o) V I m p + k V 2 <|) m ’ m f OT p 


where the coefficients y are evaluated from the recurrence relation 


(3.95) 
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Tot + 2/2 Tot -1/2 ^>m^m 


(3.96) 


In the above equation A denotes the cell face area and v denotes the cell volume 
and subscripts w, e, s, n, P denote the west, east, south, north, and cell center, 
respectively. Figure 3-1 shows a sketch of the control volume. 



The in-scattering term on the right hand side of this equation contains the phase 
function 0 which can be prescribed for both isotropic and anisotropic scattering by 
the medium. For a linear anisotropic scattering, the phase function may be written 
as. 


4Wm = 1 • 0 + [Fot IV + ^^m- +T lm T l «] 0.97) 

where a 0 is the asymmetry factor that lies between -1 and 1. The values -1,0,1 denote 
backward, isotropic and forward scattering, respectively. 
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Equations (3.95) involves 5 intensities out of which 2 intensities are known from 
the boundary conditions. The other cell face intensities can be eliminated in terms 
of the cell center intensity by a spatially weighted approximation. For example, if 
the east and north cell face intensities are unknown, they can be eliminated in terms 
of cell center intensity by using 


/ -(2-cc)/ 

V = - 1 


a 


i =— 

a 


(3.98a) 


(3.98b) 


where a is the weighting factor. a=l yields a upwind differencing scheme and oc-1/2 
yields a central differencing or the "diamond differencing" scheme. Substituting 
Equation (3.98) into the integrated discrete-ordinate equation results in a algebraic 
equation for the cell center intensity. 


^m AI m a + ^m BI m, + S l- S 2 

Imp ~ MN + Us + «( h(I ] V 


(3.99) 


where 


A = A e (1 - a) + A w a 
B = A n ( l-a) + A s a 
S 2 =a*VZ ¥ 
a 

S2 = ot V w m ' Z m 'p 


Equation (3.99) is appropriate when both the direction cosines are positive and the 
direction of integration proceeds in a direction of increasing x and y. For negative 
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direction cosine, the direction of integration is reversed and the integration sweep is 
started from the appropriate corner of the domain. Thus, for each set of direction 
cosines, the intensities at all the cell centers are obtained by marching in the 
appropriate direction. This procedure is repeated for all the ordinate directions. The 
in-scattering term is evaluated explicitely by using the previous iteration values, 
thus decoupling the m ordinate equations. 

Assuming local radiative equilibrium, the source term for each computational cell 
i.e. the net radiative heat flux, is given by. 


S = Vk2,I m w m -V4IlkI b 

m 


( 3 . 100 ) 


This source term is added to the gas phase energy equation. The gas phase equations 
and the radiative transfer equations are solved in an iterative manner. 
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4. GEOMETRY AND ITS COMPUTATIONAL REPRESENTATION 

The REFLEQS code uses a single domain structured grid approach. One important 
feature of the code is the use of a non-orthogonal body-fitted grid with internal 
blockages. As a result, many flow problems with complex geometries can be 
handled. In this section, a brief description of the geometry, the Body-Fitted 
Coordinate (BFC) system, and the internal blockage concepts are presented. 

4.1 Geometry and Grid 

The REFLEQS code can be used to solve any fluid flow problems with complex 
geometries as long as the flow domain is covered with a structured grid. A grid is a 
discrete representation of the flow geometry and computational domain. A grid is 
considered to be structured if there exists three grid lines (for a 3D problem) to 
identify three distinctive directions and any face of a control volume is on two grid 
lines. In other words, a single (I,J,K) index can be used to identify a cell or a point of 
a structured grid. Since the REFLEQS code allows internal blockages, the above 
mentioned structured grid is able to cover most flow problems despite its geometric 
limitations. An example of a two-dimensional flow problem is shown in Figure 4- 
la and two possible grid layouts, one with blockage and another without blockage, 

are displayed in Figures 4-lb and 4-lc. 

4.2 Bndy-Fitted-Cnordinate (BFC) Sy stem 

In the REFLEQS code, three different coordinate systems are used: the Cartesian grid 
system; 2D axisymmetric polar system; and an arbitrary Body-Fitted-Coordinate 
system. The Cartesian and polar systems are specialized grid systems for problems 
with simple geometries and they could offer some savings of the storage and 
execution time. In general, the REFLEQS code works with the BFC system and on a 
BFC grid. Figure 4-lb and 4-lc are just two simple examples of such a grid system. 
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Figure 4-1. 



(c) 

A Sample Two-Dimensional Flow with Two Possible Grids with a 
without Blockages 


42 




lUmilJ 


Mathematically, a BFC system can be viewed as a coordinate transformation from 
physical domain to computational domain as illustrated in Figure 4-2 for a 2D 

problem. The % and T] coordinate lines and the I and J indices have a direct 
correlation. In order to identify each point on the physical domain, a base Cartesian 
coordinate system is always used as shown in Figure 4-2a. As a result, the 
transformation can be expressed as 


x = X(%,t\, Q 

£ = £(*, y,z),T|=r|(*,y,z), C = C(*,y,z) or y=Y(6 T \,Q (41) 

z = Z (£ t\,Q 

Therefore, each point on a physical domain can be identified by a triplet (x, y, z). 
This triplet is a function of ri and £ or is associated with a particular (i, j, k) index. 
Therefore, the REFLEQS code works on a BFC grid system expressed on a base 
Cartesian coordinates. Note that a physical domain is always transformed to a 
rectangular domain as shown in Figure 4-2. 

For future reference, a brief introduction of the BFC coordinate system is given 
below. A comprehensive description of BFC coordinate system can be found in the 
work of Thompson et al, (1985). For derivation convenience, & & is used to 
replace (^, r|, Q. In a BFC coordinate system, (^, & the covariant base vectors 
are defined as 



dr _ , 

^ e / ] O 


(4.2) 
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0 


4 1 0S'6 f4-2 


X 


(a) physical domain 



(b) computational domain 

Figure 4-2. An Illustration of the Transformation from Physical to Computational 
Domain 
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where “ is the displacement vector and is equal to xi+yj + zk > and contravariant 
base vectors are defined as 


->2 e 2 xe 3 _ 2 e 3 xe 2 ^ 3 eixe 2 

£ = — — ; e - J ; 8 - J 


( 4 . 3 ) 


where / is the Jacobian defined as J = (e 2 x e 2 ) • e 3 . Basically, e / is a base vector along in- 
coordinate line and ? is a base vector normal to the surface formed by £; and 
e k (i jtj * k). But note that £ ; - and z } are not unit bases and unit bases can be defined as 


and 



( 4 . 4 ) 



(no summation on j) 


( 4 . 5 ) 


where hj is usually called scale factors. 

It is easy to show that the covariant and contravariant base vectors satisfy the basic 
relationship 


£.-2 = d ~ 

fc i c v ij 


( 4 . 6 ) 


where Sq is the Kronecker delta. 

In the BFC system, the gradient, divergence, curl, and Laplacian operators can be 
expressed in conservative form as: 


Gradient: 



Divergence : 
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Laplacian : 


v 2 /= 


i a 2 


J d?d§ 



The geometric meaning of the above quantities can be easily explained on a 2D local 
point P curvilinear coordinates shown in Figure 4-3. For a BFC grid, the Jacobian is 
simply the volume of each cell. 



Figure 4-3. Geometrical Meaning of Covariant and Contravariant Bases for a 2D 
BFC System 

4.3 Blockage Concept 

In many engineering problems the boundaries of the domain are very complicated 
and there can be internal flow obstacles. To facilitate the use and generation of a 
single domain structured grid, the REFLEQS code employs the so called ''blockage" 
concept. With blockage capability, several parts of the flow domain can be real solids 
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sr 3 

B 

m 

or imaginary dead regions which are blocked-off. Such blocked-off regions are 
illustrated in Figure 4-4 for a 2D flow where both real and imaginary blockage is 
used. As it can be seen, a structured grid is very difficult to generate without 
blockage capability and blockage offers convenience in grid generation. The blocked- 
off regions will participate in the grid generation and numerical solution process but 
will have solid wall properties. The numerical treatment of the blockages is not 
discussed here and interested readers should consult the book by Patankar (1980). 



Figure 4-4. Illustration of the Use of Blocked-Off Regions 
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5. DISCRETIZATION METHODS 


The fluid flows are governed by several physical conservation laws and these laws 
can be written as Partial Differential Equations (PDE's) as presented in Section 2. A 
numerical method to solve these PDE's consists of the discretization of the PDE s on 
a given grid, formation of corresponding linearized algebraic equations, and the 
solution of the algebraic equations. This way, a final set of discrete numbers on a 
grid is obtained which represents the numerical solution of the PDE's. In this 
section, the discretization of the governing equations is presented. The finite- 
volume approach is adopted due to its attractive capability of preserving the 
conservation property. In the following, the storage locations of the dependent 
variables are first discussed. The discretization process then follows in more detail. 

5.1 Staggered versus Colocated Grid Approach 

The staggered grid approach was widely used for incompressible flow simulations in 
the past. With the staggered grid approach, proposed by Harlow and Welch (1965), 
the velocity components are stored at positions between the pressure nodes as 
illustrated in Figure 5-la. Such an approach ensures the pressure is readily available 
for momentum equations and velocity components are available for the continuity 
equation without interpolation. As a result, a proper pressure coupling is 
guaranteed and the well-known checkerboard instability is prevented. However, 
the disadvantages of the staggered grid approach are well known and the following 
are just two examples. 

• It is not easy to extend the staggered grid approach to non-orthogonal 
curvilinear grids. Several proposed extensions are extremely complex 
to apply and can cause loss of accuracy. 

• Many state-of-the-art CFD methodologies are difficult to apply with a 
staggered grid approach such as multigrid, local grid refinement, and 
multi-zoning. 

The current state-of-the-art approach is based on the colocated grid (or non-staggered 
grid) proposed by Rhie and Chow (1983) and later by Peric (1985). The grid 
arrangement of this approach is illustrated in Figure 5-lb. This colocated grid 
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approach has many advantages over the staggered grid approach and is used in the 
REFLEQS code. This approach evaluates the cell face velocity using momentum 
equations and consequently the cell face velocity is directly linked to a third pressure 
derivative term. This linkage ensures that the checkerboard instability is 
eliminated. With a colocated grid approach and Cartesian velocity as dependent 
variables, many coding complications are avoided and more accurate solutions can 

be obtained, particularly for viscous flows. 



Figure 5-1. Illustration of Variable Storage Locations for: (a) Staggered; and 
(b) Colocated Grid 



It is noted that all governing equations, except continuity equation, can be expressed 
in a general form as 
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where <p can stand for Cartesian velocity components, total enthalpy, turbulence 
kinetic energy, mass concentration, etc. P is the effective diffusivity and is the 

source term. Therefore, the above equation, in convection-diffusion form, will be 
considered for discretization. The continuity equation will be discussed later in 
Section 5.7. 

First of all. Equation (5.1) needs to be transformed to BFC coordinates using the new 
independent variable t,(x, y, z ) , T](x, y, z) , and £(x, y, z). Without detailed derivation, 
it is sufficient to write down the transformed equation in (t, , T), Q coordinates as 
follows: 


ttw+iiW * 


y-e 




pjg 1 


ft 


d(|) 

a?j 


+ 


(5.2) 


with 


ik -*i -it 

g 1 = £ • e 


(5.3) 


The discre tiz ation involves an integration of Equation (5.2) over a control volume 
as shown in Figure 5-2. That is: 


[Equation (5.2)} J d %d rj d £ 


(5.4) 


Note that At, = At] = AC, = 1 has been chosen on every control volume with the 
REFLEQS code. This is purely for computational convenience. In the next several 
sections, individual terms of the above integration will be discussed. 
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Figure 5-2. The Labeling Scheme of a Control Volume 
5.3 Transient Term 

Consider the discretization of the transient term of Equation (5.2), i.e., 


T= (7 P$)d%dT\ d£- P< ^ (55) 


In the above, superscript "o" stands for a value at an old time, variables without "o" 
supersdpt are at the new time, and V stands for the volume. Note that the above 
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discretization holds true for both Euler first-order and Crank-Nicholson second- 
order schemes. 

5.4 Convection Term and Different Convection Schemes 


By defining a physcial contravariant velocity component jJ such that v ~ u e k , the 
convection term can be rewritten as 




1 d 


J-pU 1 * 


l_d_ 
+ Jdn 


1 pU 2 i D 


h 


l d_ 
+ JdC 


J 


pU 3 0 


(5.6) 


Integration of the above term over a control volume gives. 


C = G e - C w + C n - C s + C h - C l = G A " G w^w + G r$n ' G A + G A ' G M 


(5.7) 


with G e defined, for example, as 



(5.8) 


where G's represent the mass flux through a face of the control volume. The 
subscripts, e, w, n, s, h, and / represent the values evaluated at east, west, north, 
south, high, and low faces as shown in Figure 5-2. The evaluation of the mass 

fluxes G, will be described in Section 5.7 and the evaluation of 0 values at control 
volume faces are described next. 

Without loss of generality, the evaluation of (f> e is discussed on a uniform grid as 
shown in Figure 5-3. The following various schemes can be used for this purpose: 
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n 



Figure 5-3. An Illustration of the Cell Face Value Evaluation 
5.4.1 First-Orde r Upwind Scheme 

First-order upwind approach will evaluate fa using either fa or fa depending on the 
flow direction at point e. Mathematically, it can be expressed as 



GA ,ifG e >0 

GJ> E ifG,<° 


or in more convenient form as 




(5.9) 


(5.10) 


This scheme has first-order accuracy and is one of the most stable schemes. 

5.4.2 Central Difference Scheme 

Pure central difference approach will evaluate fa by averaging the values at point P 
and E. That is 
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(5.11) 


■CN_ r *E + $P 
e ” 2 


It is widely known that pure central difference scheme may cause unphysical 
oscillations. For most problems, some damping (or artificial viscosity) is needed for 
stability. In practice, central difference scheme with damping is constructed as 

C, = d,C* n ’ + (l-<i,)C?' (5.12) 

with d e , the damping coefficient, representing the fraction of Upwind scheme used. 
Equation (5.12) can be re-written as 


C e = G e 


4>p + <[>£ 


-d. 




(5.13) 


Comparing Equation (5.10) with Equation (5.11), it is dear why the scheme is called 
central difference plus damping. d e = 0 yields the pure central scheme, while d e = 1.0 
results in the upwind scheme. In this report, central difference scheme often refers 
to pure central plus the damping term as given above. 

5.4.3 Smart Scheme with Minmod Limiter 

It is observed that the damping coefficient is constant for a central difference 
scheme. There are many flow situations where damping is needed only in certain 
limited regions. Therefore, a smart scheme is designed to adaptively calculate 
damping coefficient depending on the local flow property. The minmod limiter can 
be used to obtain damping coefficient d e as 

d e = minimode (l, y e ) (5.14) 


with 
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(5.15) 


_ §E-U~$P-U 

Ye “ 4 > e - <1 > P 


and U = sign 



The MINMOD function is defined as: MINMOD (a, (3) = sign (a) max [o, min (j a\ 

( 3)1 

This scheme will reduce to upwind scheme if there exists local extrema such as 
across shocks (y e < 0)/ to central difference scheme when y e ^ 1, and to second-order 

upwind scheme for 0 < y e < 1- 

5.4.4 Other High -Order Schemes 

For compressible flows with shocks, several high-order schemes with limiters are 
very accurate for shock capturing. Three such schemes are introduced in the 
following. 

Osher-Chaki-avarthv Scheme: Osher-Chakravarthy scheme evaluates the damping 
coefficient as 


d e = minmod (/?, yj + minmod (l, £ y e ) (5.16) 

1 a a _ 3 ' Tl 
with R = j and p - 

Roe's Superbee Scheme: Roe's Superbee scheme is used to obtain the damping 
coefficient as 


d e = Max [o, min (2y, 1), min (y, 2)] (5.17) 

van beer's MIISCI. Scheme: van Leer's MUSCL scheme is used to obtain the 

damping coefficient as 
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5.5 


Diffusion Terms 


The diffusion term in Equation (5.2) can be split into two parts: main diffusion 

(j = k); and cross diffusion (j ^ k). Let's consider main diffusion first, i.e., 




1 a 




Tlx 1 * — 
18 *5k 


k = l,2, or 3 


(5.19) 


Without loss of generality, k = 1 term will be used for derivation. Integration o i&M 
term over a control volume, as shown in Figure 5-2, leads to 


Dhjd&ndC-- 




r Jg 


a 


5(j) 

K 


(5.20) 


It is easy to show that 


t li _ A? _ A 

J ~ h} sin cq (5.21) 


where A is the area of the control volume face along the % direction and ai 

— » — * 

represents the angle between vector e j and the plane formed by Qj srid £y By 
defining 


r -a 

De — -j . 

$ h 2 sin ctj 


(5.22) 
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we have 


JJJ Di,;iM=D' 5 (t E ^)-Dj(* r * w ) = -(D i {t Dj)t, t D^ t + Df+w (5.23) 

Therefore, D^, the main diffusion coefficient needs to be evaluated at the faces of 
each control volume. 

The cross diffusion term is written as 


& m l± 

c JHk 


IT g‘ 


& 


3d 

a^J 


, j 


(5.24) 


Consider D* 1 for illustration purpose: 


f f ^ (5.25) 

mj J J 

By defining D n =jlg n and assuming “ 4 f^P + + + ^ne). etc., we have 

JJJ Df / Df 1 (<t> N + <t>NE - 4>s - fe) - Di, 1 + <>NtV - <)>S - <l>sw) ®. 26 ) 

Other cross terms can be similarly discretized. 

5.6 Final Linear Equation Set 

From Subsections 5.3 to 5.5, the general convection-diffusion Equation (5.2) has 
been discretized term by term on a control volume. Assembling all the terms 
together, the final linear equation is resulted and can be written in general form as: 
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A$p = A + A$ e + A^ s + Aj$ n + A l $ l + Aj$ h 

+ ^ sw^sw + a se$se + a nv$nw + a n$ne 

+ A L$LS + A LN$LN + A HS$HS + A Hh$HN 
+ A WL$WL + A WE$WH + A E L$EL + A Eh$EH 

+ S* 


5.7 Discretization of Mass Cons e rvation and Mass Flu* Evaluatio n 

As mentioned before, mass conservation equation is a special one which can not be 
written as the general convection-diffusion form (Equation (5.2)). Moreover, in a 
pressure-based method mass conservation is used to determine pressure field. For 
this purpose, the mass conservation equation needs special attention. The general 
mass conservation equation can be transformed into BFC coordinates as 

7 d 4 + mS ,p ^ =0 (528) 

Therefore, integration over a control volume yields 


pV - p°V° 


Af 


+ G - G,„ + G„ - G c + Gl 


w 


G[ = 0 


(5.29) 


Wjtfi ^ e \hi ^ ^ j g for example. Note that G's are the mass flux through a control 
volume face as mentioned in Section 5.4. 

Our next task is to evaluate mass flux G's at control volume faces such that the 
checkerboard instability is eliminated for a colocated grid arrangement. This is 
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jan 


achieved by averaging momentum equation to the cell faces and relating the cell 
face velocity directly to the local pressure gradient. This procedure is presented next. 

The discretized x-momentum equation can be expressed at cell P as 


% n 


a 


u 

V 


+ 



nb 




u°„ + S“ 


( 5 . 30 ) 


where nb refers to all the neighboring cells and the transient term is explicitly 
expressed to ensure that the momentum equation at cell face includes the effect of a 

time step. Dividing the above equations by a“ yields 


where 


nb v 



+cd pK 



( 5 . 31 ) 


c = 


_P 

At ' 



( 5 . 32 ) 


The above momentum equation is at P-cell but in reality it applies to every point. 
Therefore, at cell face f, for example, we have 



v u nb 

L u 
\nb a p 


U 


nb 




if 



+ cd'fu c f + 


f cu\ 


\ a vJf 


( 5 . 33 ) 


Since we do not know how to evaluate for example, the following quantities at 
/ will be obtained by averaging the same quantity from two neighboring nodal 
points. 
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(5.34) 


“f- 


( ^ 

7 nb 


v-> -no 

^~7 U nb - 


(s u « 


V”‘ “p if \ r If 


Therefore, by using Equation (5.31) we have 


U r u 

a nb S u 

L — u nb + — 

\_nb a p 


Pi 


1/ 


„u C u 

^nb 

Z— U nb + — 

\nb a p a P J 




(5.35) 


where overbar means average from point P and E. Hence Equation (5.33) can be 
written as 


(l + cdf) u , = {2 + cd“) + rf“ (g) - cd“ - d; (!) 


+ CdJ u° f 


(5.36) 


and the above equation is further reduced to 


Uf=u p + d“ 


W\ ffy 


telp [dx/fj 


■ ci “ (“/-"?) 


(5.37) 


with 


d u 



1-cd 


u 

V 


(5.38) 


Now it is clear from the above equations that the cell face velocity is obtained from 
an average of the two neighboring point velocity plus a second-order pressure 
correction and a second-order previous time velocity correction. The pressure term 
serves as the mechanism to avoid the checkerboard problem and the previous time 
velocity serves as a mechanism to obtain a time-independent steady solution. 
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The y-component and z-component velocities can be similarly evaluated 
contravariant velocity component at east face can be thus evaluated as. 

Uj=u r F lx + v r F ly + w r F lz 

with 

F lx = h l^-^' F ly = h l el -1 ; F lz = h l el -k 


and the 


( 5 . 39 ) 


( 5 . 40 ) 
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6. VELOCITY PRESSURE COUPLING 


In Section 5, the mass conservation equation and other governing equations are 
discretized on an arbitrary control volume. It is noted that both pressure and 
Cartesian velocity components appear in the mass conservation equation and the 
main variable is not clearly identifiable. It is the task of this section to reformulate 
the mass conservation equation such that the pressure increment (or pressure 
correction) is the main variable. As a result, the velocity and pressure fields are 
linked together. This approach is classified as a pressure-based or pressure 
correction method. 

For ease of presentation, only incompressible flow and Cartesian grids are used. 
Extension to BFC and compressible flows are straightforward. The discretized 
equations for mass and momentum conservation can be re-written as 


'Luf-O 

f 



u i = 'LA i nb u irnb -D i V i p + S 

nb 


( 6 . 1 ) 


y Y 

where " stands for the summation over all faces, 

f no 

neighboring cells. 


is the summation over all the 


In the following, three velocity-pressure coupling algorithms are described which 
are available from the REFLEQS code. 


6.1 SIMPLE Algorithm 

The SIMPLE algorithms (Patankar et al , 1972) relies on using a guessed value for 
pressure p*, and using the momentum equation can be used to obtain the velocity. 
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u\. If the pressure has an increment p', the velocity components will also have a 
corresponding increment u\. This increment can be written as 


In the SIMPLE algorithm u[ nl are neglected and the above increment equation is 
approximated as 


A 1 u\ = -E? V/ or w, = - 



( 6 . 3 ) 


Therefore, at cell faces, the increment is 



EL 

UJ 


V /P 


We would like to enforce the new velocity, 

u f~ u *f + u f 

to satisfy the mass conservation, i,B., 


or 


-£u' = !u f 

f f 





f 


( 6 . 4 ) 


( 6 . 5 ) 


( 6 . 6 ) 


( 6 . 7 ) 
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Note that the right hand side of the above equation is simply the mass imbalance. 
The above pressure correction equation is a Poisson equation and any linear 
equation solver can be used to solve it. After p' is obtained, a new pressure and a 
new mass conserved face velocity Uf are obtained using Equations (6.4) and (6.5). 

6.2 SIMPLEC Algorithm 

SIMPLEC algorithm (van Doormal and Raithby, 1984) differs from SIMPLE in the 
approximation of the momentum increment equation. With SIMPLEC algorithm, 
it is assumed that 


^ A nb u i,nb ^ ^nb u i 
nb nb 

As a result, the increment momentum Equation (6.2) becomes. 


A;-z4>;=-D f v,. P 

i nb 


(6.8) 


(6.9) 


Since 7i is always larger than 



nb ' 


the coefficient 



nb 


( 6 . 10 ) 


is always positive. 

With Ap replacing A' p the same derivation performed for SIMPLE can be used and 
the pressure correction equation is also given by Equation (6.7). It is obvious that the 
modification from SIMPLE to SIMPLEC is very small. However, in most problems, 
the efficiency improvement is significant. 
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6.3 PTSO Algorithm 


PISO is a Predictor-Corrector algorithm and is more mathematically sound 
compared with SIMPLE and SIMPLEC since no assumptions are made to the 
momentum increment equation. The algorithm is described as follows, with a 

guessed pressure, p\ the momentum equation is first solved to get u]. This is the 
predictor step and can be written as 


Predicton 




I,A i nb ul„ b -DVi + 

nb 


s. 


( 6 . 11 ) 


Next, two correctors are performed. First-corrector will provide a new u i and p 
such that 


lu/=0 

/ 

%u? = I,A&w-lWf~+S i 

nb 


(6.12) 


and second-corrector will get the velocity w/ and pressure, p, such that 


Z«/ = 0 
/ 


(6.13) 


The similarity between PISO and SIMPLE is obvious if the new momentum 
equation is subtracted from the old momentum equation. That is: 
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First corrector: 


( 6 . 14 ) 


'+ ** 


u i = u i ~ u i 


A 1 U: =-DV-p with < ,» 

P =P -P 


** ' ** 


Second corrector: ^nb u i,nb'^^ P with u,- u i -u i ,p p p 


( 6 . 15 ) 


Since U inb is known at second-corrector, the above increment equation has the 

similar form of SIMPLE and SIMPLEC. However, it is clear that PISO does not make 
any approximation regarding the momentum increment equation. Also it is worth 
mentioning that the number of correctors does not have to be two and more 
correctors can be performed. 
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7. BOUNDARY CONDITIONS 


In order to completely specify the mathematical problem, it is necessary to supply 
the conditions at the boundaries of the solution domain for all the dependent 
variables. These conditions are usually either specification of the value of the 
dependent variable, or the value of the associated flux, or a combination between 

the two. 


REFLEQS offers a variety of boundary conditions which include: 


a. Inlet : where all the dependent variables are fixed except pressure; 

b. Exit : where all the dependent variables are extrapolated from the 

interior except pressure. Pressure can be fixed (FIXP), extrapolated 
(EXTRAP), or both (FIXPEX); 

c. Symmetry : where the normal gradient of all the variables are equal to 
zero except the velocity component normal to the boundary which is 

set to zero; 

d. Solid Wall : where the velocity is fixed, pressure is extrapolated, and 
the value or the flux is known for other dependent variables; and 

l Periodic Boundary Condition: where the two boundaries exchange the 

information. 


For each dependent variable, at all boundaries, appropriate modifications of the 
discretized equation at the "near boundary" nodes is required. For most of the 
dependent variables, the boundary conditions are implemented by modifying the 
source terms S 17^ and SPf 


u +stJ » 

nb 

' F ~ XA nb -SP^ (7.1) 

nb 


A brief introduction of the above boundary conditions are provided below. 
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7.1 


Tnlpf with Specified M ass Flow Rate (FIXMl 


At an inlet, the amount of incoming mass flow rate m JN through the boundary cell 
face and the incoming fo property should be specified. In this case, the finite- 
difference coefficient connecting the boundary node to its neighboring internal node 
is set to zero and then the SU<p and SP$ source terms are modified as: 


SUq = SUq + m IN <t> B 
~ ' m IN 

The pressure correction equation is not modified in this manner. A link coefficient 
with the boundary node is set to zero for p' -equations. 

7.2 F.xit Boundary 

At the exit, the fluid flows out of the calculation domain, and the information of 
most dependent variables is not available a priori. Fortunately, the information 
there is often not important if the convection is the dominant activity. Therefore, 
all the dependent variables except pressure can be extrapolated. This is equivalent to 
setting the boundary link coefficients to zero. As for pressure, it is fixed at the exit 
for most incompressible or subsonic Mach number flows (this is called FIXP 
boundary condition). However, for supersonic flows, pressure should be 
extrapolated at the exit (EXTRAP boundary condition). REFLEQS offers a third 
boundary condition, FIXPEX, to take care of the situation where the exit flow 
condition supersonic or subsonic is not known a priori. 


7.3 Symmetry Condit ion 

At symmetry boundary, a zero mass flux through the boundary is assumed. At the 
same time, the normal derivative of all the scalar variables is set to zero. 
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7.4 Solid Wall Boundary 


For laminar flows or turbulent flows without using wall functions, the no-slip 
assumption is made at the solid walls. As a result, the mass flux through a wall is 
set to zero. On the other hand, the pressure at the wall is not known, a priori, and 
extrapolation from the internal cells is used. For other dependent variables, a 
constant value, a constant flux, or a combination of the two is needed at the wall. 
For example, for the energy equation, either constant temperature or adiabatic 
condition can be assumed. For turbulent flows with wall function boundary 
condition, a special procedure is carried out. 


The wall functions described by Launder and Spalding (1974), are derived from 
experimental and analytical knowledge of the one-dimensional Couette flow which 
exists near the wall. A semi-empirical universal function of non-dimensional 
distance normal to the wall, y + , is: 



pfy-u T 

M- 


(7.4) 


In the above definition, Sy is the distance normal to the wall and u T is the "friction 
velocity" given by: 


«r = 



(7.5) 


In the internal sublayer (y* > 11.63) the velocity variation may be described by a 
logarithmic relationship i.e .: 


u = 


K 


tn (Ey + ) 


(7.6) 


where E = 9.70and k = 0.4034 are experimentally determined constants. 
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In both the viscous (y + < 11.63) and internal (y + > 11.63) sublayers, the shear stress is 
calculated from the product of effective viscosity j u e ff arid normal velocity gradient 
dii/dy, i.e.: 

< 7 - 7) 

where 


li for y + <11 .63 

/ + «« (7 ' 8) 
vturb ior y >n - 63 

Near the wall, the transport equation for the turbulent kinetic energy, k, reduces to a 
balance between the local production and dissipation of k to give: 



m 



= pe 


(7.9) 


The velocity gradient may be replaced from Equation (7.7) and the dissipation rate 
from: 


„ k 2 

to give: 

* w = Cf/ 2 pk = pu 2 T 

Hence, it follows from Equation (7.5) 


(7.10) 


(7.11) 
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(7.12) 


pcj[ 4 k 1/z u 

tw= L M E y + ) 


7.5 Periodic Boundary 

The periodic (cyclic) boundary conditions appear in the circumferential direction if 
the two ends of the calculation domain in the z-direction join up with one another. 
This can occur in a polar-coordinate direction in which the whole angular extent 
from 0 to 360° is to be considered (Figure 7-la) or when "repetition" is present in the 
flow pattern in the angular coordinate direction (Figure 7-lb). 

The general rule is that whenever identical conditions are to be expected at z = o and 
z = last z, the finite flux, equal on both sides, is to be expected at that surface, then the 
boundaries are cyclic. 

In this circumstance, the boundary conditions can be specified as follows: 

§LB = ' §HB = $2 ' an< ^ = ^ HB (7-13) 

where indices HB and LB denote high boundary (k = N) and low boundary (k = 1), O 
is the convective flux in the z-direction. 
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8. SOLUTION METHOD 



In the previous sections, the discretization method, the derivation of the pressure- 
correction equation, and the implementation of the boundary conditions have been 
described. The task of this section is to assemble all the elements together and 
describe how the discretized governing equations are solved. 

8.1 Overall Solution Process 

The Navier-Stokes equations are a set of highly non-linear and coupled differential 
equations. In the discretization of the governing equations, it is seen that the 
coupling among different dependent variables is one iteration delayed and the non- 
linear terms are linearized. Therefore the above approach is usually called an 
iterative method. Due to the de-coupling procedure, the final linear equation set for 
each governing equation is well defined and can be solved separately using any 
linear equation solver. This approach is called segregated or equation-by-equation 
approach. In the following, the overall solution process of the segregated method is 
described to solve transient Navier-Stokes equations. 

Prescribe an initial fluid flow field {<p °] at the beginning of the 
computation (t°). 

Make one increment in time (t n+1 = t n + At) and we try to get solution 
{Qn+i}. For this purpose, an iterative method is used. Therefore, set 

(<P n+1 k= {<pn}. 

With known (tp n+ ^)k/ all the link coefficients A n b are evaluated. 

The momentum equations are solved first to obtain {' u i } k+ j using 

the known pressure {p n+1 } k . This velocity field will not satisfy the 
mass conservation. 
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e. Pressure correction is solved to obtain the new pressure {ip /^ + j ^rid 

new velocities are obtained as K +J }* +1 . This new velocity field may 
satisfy the continuity equation but not satisfy momentum equation. 

f. Any other scalar equations, such as energy, turbulence, combustion 

equations are solved to get 

g. Steps c-f are repeated until all equations are satisfied or a 
predetermined convergence criterion is met. At this point, the 
solution at n+1 time step is calculated. 

h. Steps b-g can be used to obtain a transient solution at any time 
locations. 

8.2 Linear Equation Solvers 

In the overal solution process, one important component is to solve the linear 
equation set resulted from each governing equations. This linear equation solver 
usually consumes a large fraction of cpu time and may require large amounts of 
memory. Therefore, the selection of Unear equation solvers is always a trade-off 
among execution time, storage and accuracy needed. A direct solver such as 
Guassian Elimination is an accurate method but requires prohibitive storage and 
execution time. It is rarely a useful solver for segregated iterative approach. The 
best candidate for solving sparse banded matrix Uke ours is the iteration based 
solver. In the 60's and 70's Alternating Direction Implicit (ADI) method was very 
popular. Despite the fact that it requires very small storage and cpu time per 
iteration, ADI method becomes prohibitively slow to solve the large system. 
Therefore, ADI is also not a good candidate. In REFLEQS, two advanced solvers are 
developed. One is the Whole Field Solver (WFS) which is an enhanced Stone's 
Strongly Implicit Procedure (SIP). The WFS, unlike SIP, is free from any adjustable 
a- parameter and is symmetric. Another solver is the Conjugate Gradient Squared 
(CGS) solver which may provide an accurate solution for difficult linear system 
with reasonable storage and execution time cost. The two solvers are described in 

the following sections. 
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ri 



8.2.1 Whole F tpld Solver 

Figure 8-1 illustrates the backsubstitution node arrangement and the nomenclature 
for the WFS discussion. 



With the i,j,k indices omitted the finite-difference (FD) equation can be expressed as 


<y(> = a$ M + «w4>f_ i + a bfij+i + a ^j-i + a tftk+i + a i$k-i + su (8,1) 

The difference between WFS and SIP solvers is the solution algorithm of the z- 
symmetric 3D WFS for the system of the FD equations can be summarized as 
follows. 

a. In SIP only 3 points are treated implicitly in backsubstitution (P, H, N) 
while in WFS (P, L, H, N). 

b. Nonsymmetric arrangement of cancellation nodes is used in SIP while 
in WFS (SH, SN) symmetric arrangement is retained, (SL, SH). 

c. No extra memory storage is needed in WFS compared to SIP. 
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8 . 2.2 


Solvgi 


Conjugate gradient type solvers have become very popular in recent years since they 
often provide a faster convergence rate. They also have many other advantages 
(vectorization, no user specified parameters, etc.) over classic iterative methods. 
Among various nonsymmetric and nonpositive definite conjugate gradient solvers, 
the CGS algorithm has been found to yield faster convergence and is simpler to 
implement in the CFD code. The CGS algorithm applied to the system AX = b is 
expressed as follows: 


Initialization (n=0), 


r 0 = b - AX 0 

q 0 = P-i = 0 ; P-i = 1 


~ T o 

Pn =r o r «Pn = J- 
rn-l 

u n = r n + p n q n 

Pn = U n + Pn {tfrt + Pn Pn-l ) 
Vn =A Pn 


Interaction ( n>0 ), q n+ \ - u n ' a n v n 

r n+l 

x„ + i=x n + cc n [u n + q n+l ) 


The convergence rate of conjugate gradient algorithms depends on the spectral 
radius of the coefficient matrix and can be effectively accelerated by preconditioning 
the system. To a matrix, the linear system AX = b is transformed into an equivalent 
one whose coefficient matrix has a lower spectral condition number. This is 
obtained by using an approximate matrix p, and this transformed system is solvable 
by a direct method, and applying the CGS algorithm to the equivalent system 
AX=p A AX = p' 1 b = b- In Practice, each time in a CGS iteration one has to compute a 
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matrix-vector product of the form y = AX> with A = p A A • The vector y is actually 
computed by solving a linear system of the form py = AX. Incomplete Cholesky 
decomposition is used as a preconditioner in the code. It has been found that 
preconditioned CGS is very efficient and robust. 
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9. CONCLUSIONS 


As a result of this study, an advanced CFD- combustion code, REFLEQS, has been 
developed for analyzing two-phase spray flows in rocket combustors. State-of-the- 
art numerical schemes and solution algorithms were incorporated into the code for 
accurate flow simulations in complex geometries. Colocated grid arrangement, 
strong conservation of momentum equations and high order numerical schemes 
were used to ensure good solution accuracy. 

To simulate flow and combustion in liquid rocket engines, several physical models 
were incorporated into the code including: 

• turbulence models (Baldwin-Lomax, standard k-e, extended k-e, low 
Reynolds number, and multiple-scale models); 

• H2/O2 and hydrocarbon/ O 2 combustion models (equilibrium, and one- 

step and reduced mechanism models), 

• PDF based turbulence combustion interaction models; 

• two-phase Eulerian-Lagrangian spray dynamics model; and 

• discrete ordinate thermal radiation model. 

The code was systematically and extensively verified against a large number of test 
problems to ensure solution invariance to coordinate orientation, grid skewness, 
boundary condition order, initial conditions, etc. Results of the testing, validation, 
and application studies are presented in the REFLEQS: Valiation and Applications 
Manual The code has also been successfully applied for several liquid rocket engine 
configurations including SSME, AMCC, STME, and OMV. It has also been used for 
several turbomachinery applications, including NASA MSFC's Bearing Seals and 
Materials Tester (BSMT) flow problems. 
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10. RECOMMENDATIONS 


Based on the experience of work during this project and assessment of current 
trends and future needs, several improvements are recommended for REFLEQS. 
The suggested enhancements are classified in the following three categories: 

a. Physical Model Enhancement; 

b. Numerical Method Improvement; and 

c. Code Validation and Applications. 

10.1 Physical Mo del Enhancement 

To resolve complex spray combustion dynamics in rocket engines, the following 
physical models need further enhancement: 

Multi-Step, Finite-Rate, Chemical Kinetics Model: The chemical equilibrium 

analysis has been widely used to predict species and temperature distributions in the 
combustion chamber and the nozzle. However, this analysis is valid only under 
low Mach number and high temperature conditions for which the chemical time 
scale is much smaller than the flow time scale. In regions of low temperature (near 
cooled walls or regions where the mixture ratio is far from stoichiometric) or high 
speed flow (i.e., in the nozzle region), the equilibrium assumption is no longer valid 
and the effect of finite rate chemistry must be taken into account. The model should 
be formulated in the general form and demonstrated for established reaction 
mechanisms (such as 7-step O 2 -H 2 and multi-step C n H m - O 2 ). 


Supercritical Spray Model: The properties of fluids change drastically near the 

critical point. Most rocket engine thrust chambers operate under supercritical 
conditions. Advanced models, such as those being developed at Pennsylvania State 
University and University of Tennessee Space Institute (UTSI), should be 
incorporated to include the effects of supercritical phenomena on the vaporization 
process under these conditions. 

Coupled Atomization and Spray Dynamics Model: Existing spray dynamics models 
require droplet size distribution and spatial resolution for droplet formation rate. 
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This information, assumed a priori, has great influence on spray dynamics results. 
Recently, advanced atomization models, based on Jet Embedding and VOF methods, 
have been developed (Giridharan et al, 1991 and Giridharan et al, 1992). The most 
suitable models for analysis of single injector and multi-injector interactions should 
be incorporated and numerically optimized for shear coaxial, shear/swirl, and 
impinging injectors. 

Conjugate Heat Transfer Model: A conjugate heat transfer model should be 

implemented to assess the effects of heat transfer between various engine 
components (such as the heat flux across the injector face to the incoming cryogenic 
propellants). The conjugate heat transfer model needs to be incorporated in 
conjunction with the multi-domain capability. 

10.2 Numerical Algorithm Im provements 

Current trends in design methodology indicate a shift from analyzing component 
performance to a system of component interactions. Multi-domain, multi-media 
simulation with a large number of cells (106) are needed. To achieve this, the 
following enhancements are recommended for incorporation into the code. 

Multi-Domain Capability: The analysis of interior and exterior flow in rocket 

engines can be greatly optimized by the use of multi-domain capability. It allows 
efficient grid generation and limits the number of blocked computational cells in 
the domain. It is also useful for grid refinement in localized regions of high 
gradients. A general purpose multi-domain capability should be incorporated into 
the code. 

Adaptive Gridding Capability: This feature allows the grid to adapt itself to the 
solution of the problem. Therefore, the grid points are automatically clustered in 
regions of high flow gradients and become relatively coarser in regions with small 
gradients. This facilitates the efficient utilization of the computational cells in the 
flow domain. A moving grid capability should be used to implement the adaptive 

gridding feature into the code. 
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Convergence Enhancement Capability: Advanced methods such as multi-grid and 
coupled-variable conjugate gradient solvers should be implemented to improve the 
computational efficiency of solving large systems of linear and nonlinear equations. 
In addition, a second-order convergence rate scheme such as the Newton iteration 
method should be implemented and assessed. These features are especially useful 
to accelerate convergence of fine grid solutions. 

10.3 C nde Validation and Applications 

The REFLEQS code can be applied to a wide variety of, subsonic, and supersonic 
flows, with and without heat transfer and combustion. 

Even though the code has been extensively checked out for a large number of 
problems (as described in the REFLEQS: Validation and Applications Manual), 

additional detailed validation studies need to be performed for cases where reliable 
experimental data are available. More emphasis should be placed on 3D validation 
efforts and on flows coupled with complex physics (spray, combustion, radiation, 
turbulence, etc.). 
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